subroutine potcalc(x,V,nx,it,F,radius,Ef)
! double precision function V(angle)
! Função para calcular o potencial unidmensional de um sistema
! Pode ser modificada de acordo com o sistema
! Requer apenas a entrada da posição angle e retorna o valor de V(angle)

! *** Sem Variáveis Implícitas
	implicit none
! *** Declaração de Constantes e Variáveis

! * Constantes
! 'nw' é o número de poços
    integer, parameter :: nw = 1
! constantes passadas pelo programa que chama essa subrotina
	integer nx
	double precision dx, Ry, a0, am, am1, am2, ec, pi, hc
! * Variáveis
! 'V' é o potencial total
! 'Vstatic' é o potencial devido ao poço quadrado
! 'Vf' é o potencial devido ao campo elétrico variando com o tempo
	integer i, it, iw
	double precision, dimension(nx) :: x, V, Vangle, Vstatic, Vf, Vfangle
	double precision V0, radius, width, dist, F, convfactor, tau, t0, omega, w, erfc, Ef
	double complex :: img, dt, cerf, cerfc, CERARG1, CERARG2, CERARG3, CERARG4, CERRES1, CERRES2, CERRES3, CERRES4, VFEXP, VF1, VF2, VF3, VF4
	real*8 rdt

common / constants / dx, dt, img, Ry, a0, am, am1, am2, ec, pi

        open(unit=10,file='anglewi.dat')
        open(unit=11,file='anglewf.dat')

! Distância entre os poços
dist = 100.D0/a0
V0 = 262.D0/Ry
!rdt = DREAL(dt)
hc = 1.054589D3/1.602189D0

! Valores a serem usados no Potencial do Campo Elétrico
t0 = 1000.D0 ! tempo no qual o valor do potencial devido ao campo é máximo (em fs)
tau = 100.D0 ! largura da gaussiana no tempo (em fs)
convfactor = 1.D-2 ! fator para converter a unidade do campo de kV/cm para mV/Angstrom
! a frequência é determinada pela razão (Energia do fóton)/hc
!w = 1.85487D-2  ! frequencia angular da onda de luz
w = (Ef)/hc

! Cálculo do Potencial -- Poço Quântico
do i=1,nx
        Vstatic(i) = F * convfactor * radius * cos(x(i))
!        write(*,524) x(i)/pi, Vstatic(i)*Ry
enddo

do i=1,nx
!        write(*,526) x(i)/pi, F, radius, cos(x(i)), F*radius*cos(x(i))*Ry
enddo

! Cálculo do Potencial -- Campo Elétrico variando com o tempo
do i=1,nx
!	CERARG1 = ((img*omega*(tau**2))+t0-(it*cdabs(dt)))/(sqrt(2.D0)*tau)
!	CALL CERROR(CERARG1,CERRES1)
!	CERARG2 = ((img*omega*(tau**2))-t0+(it*cdabs(dt)))/(sqrt(2.D0)*tau)
!	CALL CERROR(CERARG2,CERRES2)
!	CERARG3 = ((img*omega*(tau**2))+t0-(it*cdabs(dt))-dabs(rdt))/(sqrt(2.D0)*tau)
!	CALL CERROR(CERARG3,CERRES3)
!	CERARG4 = ((img*omega*(tau**2))-t0+(it*cdabs(dt))+dabs(rdt))/(sqrt(2.D0)*tau)
!	CALL CERROR(CERARG4,CERRES4)
!	VFEXP = (EXP(-((omega**2)*(tau**2))/(2.D0)))
!	VF1 = ((img*sin(omega*t0)) + cos(omega*t0)) * CERRES1
!	VF2 = ((img*sin(omega*t0)) - cos(omega*t0)) * CERRES2
!	VF3 = ((img*sin(omega*t0)) + cos(omega*t0)) * CERRES3
!	VF4 = ((img*sin(omega*t0)) - cos(omega*t0)) * CERRES4
!	VFEXP =
!	VF1 =
!	VF2 =
!	VF3 =
!	VF4 =
!	if(isnan(DREAL(VFEXP)))then
!		write(*,*) 'VFEXP estourou'
!	endif
!	if(isnan(DREAL(VF1)))then
!		write(*,*) 'VF1 estourou'
!	endif
!	if(isnan(DREAL(VF2)))then
!		write(*,*) 'VF2 estourou'
!	endif
!	if(isnan(DREAL(VF3)))then
!		write(*,*) 'VF3 estourou'
!	endif
!	if(isnan(DREAL(VF4)))then
!		write(*,*) 'VF4 estourou'
!	endif/
!	Vf(i) = F * convfactor * x(i) * (a0/Ry) * (tau/cdabs(dt)) * ((SQRT(2.D0*pi)/4.D0)*(EXP(-(omega**2)*(tau**2)))*DREAL( (EXP(-img*omega*t0)))*(EXP(2.D0*img*omega*t0)*CERROR(((img*omega*(tau**2))+t0-it*cdabs(dt))/(SQRT(2.D0)*tau)) - CERROR(((img*omega*(tau**2))-t0+it*cdabs(dt))/(SQRT(2.D0)*tau)) - EXP(2.D0*img*omega*t0)*CERROR(((img*omega*(tau**2))+t0-it*cdabs(dt)-cdabs(dt))/(SQRT(2.D0)*tau)) + CERROR(((img*omega*(tau**2))-t0+it*cdabs(dt)+cdabs(dt))/(SQRT(2.D0)*tau))) )
!	Vf(i) = F * convfactor * x(i) * (a0/Ry) * (tau/cdabs(dt)) * (SQRT(2.D0*pi)/4.D0) * DREAL (VFEXP * ((VF1+VF2)-(VF3+VF4)))
	Vf(i) = F * convfactor * radius * x(i) * (a0/Ry) * (tau/cdabs(dt)) * 0.5D0 * EXP(-(0.5D0)*w*((2.D0*img*t0) + ((tau**2)*w))) * SQRT(pi/2.D0) * (-1.D0 + CDEXP(2.D0*img*t0*w) * (-CERF(((it*cdabs(dt)) - t0 - (img*(tau**2)*w))/(SQRT(2.D0)*tau)) + CERF((cdabs(dt) + (it*cdabs(dt)) - t0 - (img*(tau**2)*w))/(SQRT(2.D0)*tau))) + CERF((cdabs(dt) + (it*cdabs(dt)) - t0 + (img*(tau**2)*w))/(SQRT(2.D0)*tau)) + CERFC(((it*cdabs(dt)) - t0 + (img*(tau**2)*w))/(SQRT(2.D0)*tau)))
!	Vfx(i) = F * convfactor * x(i) * (a0/Ry) * SQRT(pi/2.D0) * (tau/cdabs(dt)) * (-ERF((it*cdabs(dt)-t0)/(SQRT(2.D0)*tau)) + ERF((it*cdabs(dt)+cdabs(dt)-t0)/(SQRT(2.D0)*tau)))
enddo

do i=1,nx
!        write(*,526) x(i)/pi, F, radius, cos(x(i)), Vf(i)
enddo

do i=1,nx
	V(i) = Vstatic(i) + Vf(i)
enddo

522 FORMAT(2F16.3,10F16.9)
523 FORMAT(4F16.3)
524 FORMAT(F16.3,E16.6)
525 FORMAT(3F16.3,10F16.9)
526 FORMAT(4F16.6,2E16.6)
return
end subroutine potcalc
